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Chapter  1 


Code  Validation:  Numerical  Simulation  of  Flow 
Transition  in  a  Mach  4.5  Flat-Plate  Boundary  Layer 

1.1  Introduction 

Recent  advances  in  supersonic  and  hypersonic  aerospace  technology  have  revived  the  interest  in 
boundary  layer  transition  for  high  Mach  number  flows.  A  thorough  understanding,  prediction,  and 
control  of  boundary  layer  transition  at  high  Mach  number  is  crucial  in  the  design  of  aerodynamic 
vehicles.  For  example,  skin  friction  and  aerodynamic  heating  are  considerably  lower  for  laminar 
flows  at  high  speed  than  turbulent  flow.  However,  the  mechanism  of  transition  in  high  speed 
boundary  layer  is  still  poorly  understood.  On  the  contrary  to  a  relatively  comprehensive  picture  of 
transition  scenarios  in  incompressible  flows,  nonlinear  eifects  responsible  for  transition  at  high  speed 
are  still  very  much  a  mystery  (Erlebacher  &  Hussaini,  1990).  Owing  to  the  enormous  difficulties 
of  controlled  experiments  for  high-supersonic  flows,  experimental  results  revealing  detailed  flow 
phenomena  at  the  nonlinear  stage  of  supersonic  transition  are  hardly  found  in  literature.  Instead, 
they  have  been  mostly  restricted  to  the  domain  of  linear  disturbance  growth.  Because  of  lack 
of  experimental  data  there  is  a  strong  interest  in  numerical  simulation  techniques  for  studying 
transition,  direct  numerical  simulation  (DNS)  and  large  eddy  simulation  (LES)  are  now  major 
tools  to  play  an  important  role  on  providing  understandings  to  high-speed  flow  transition. 

The  investigation  of  supersonic  transition  by  experimental  measurement  has  been  proved  to 
be  very  difficult.  Uncertainty  of  data  obtained  using  hot-wire  anemometry  is  considerably  higher 
than  that  obtained  in  subsonic  flow.  It  is  also  very  difficult  to  obtain  accurate  measurement  near 
the  wall  due  to  lack  of  spatial  resolution.  Alternatively,  the  investigation  can  be  launched  through 
theoretical  and  computational  study. 

The  theoretical  and  computational  tools  for  analysis  of  stability  and  transition  include:  com¬ 
pressible  linear  stability  theory  (LST){Mack,  1984);  secondary  instability  theory  (SIT)  (Herbert, 
1988;  Ng  &  Erlebacher,  1991);  parabolized  stability  equation  (PSE)  methods  (Herbert  &;  Bertolotti, 
1987;  Chang  &  Malik,  1993);  temporal  direction  numerical  simulation  (DNS)  (Erlebacher  Sz  Hus- 
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saini,  1990;  Pruett  &  Zang,  Adams  &  Kleiser,  1996);  spatial  DNS  (Thumm  et  a/.,  1990;  Pruett 
et  al.^  1995a,  1995b);  temporal  large  eddy  simulation  (LES)  (El-Hady  &  Zang,  1995);  spatial  LES 
(Ducros  et  al,  1996); 

In  the  linear  stability  theory,  the  base  flow  is  assumed  to  be  parallel  in  the  streamwise  direction. 
The  Navier-Stokes  equations  are  linearized  by  formulating  the  disturbance  quantities  in  the  normal 
mode  form.  This  leads  to  a  set  of  ordinary  diflferential  equations.  The  system  has  been  transformed 
to  an  eigenvalue  problem.  The  first-mode  instability  is  analogour  to  the  T-S  instability  in  incom¬ 
pressible  flow.  At  Ma  >  3.8,  the  first-mode  instability  is  purely  inviscid (Pruett  et  al^  1991).  Unlike 
T-S  instability,  the  first-mode  instability  in  supersonic  flow  are  most  unstable  when  oblique,  while 
in  incompressible  flow,  according  to  Squire’s  theorem  the  two-dimensional  mode  is  more  unstable 
than  the  three  dimensional  mode. 

In  the  secondary  instability  theory,  the  base  flow  is  defined  as  the  sum  of  mean  flow  and  the 
superimposed  primary  instability  wave.  Onset  of  secondary  instability  is  quite  sensitive  to  the 
amplitude  of  the  imposed  primary  disturbance.  There  are  three  types  of  secondary  instability: 
fundamental  (K-type),  detuned,  and  subhamonic  (H-type).  At  large  primary  amplitude,  the  fun¬ 
damental  secondary  instability  is  more  favored  (Pruett  et  al^  1991),  where  the  A-shaped  vortices  are 
aligned  along  their  peaks  in  the  streamwise  direction,  repeating  every  primary  wavelength.  Small 
to  moderate  amplitudes  of  the  primary  disturbance  often  bring  out  the  subharmonic  secondary 
instability  (Pruett  et  al^  1991),  where  the  A-shaped  vortices  are  staggered  in  the  streamwise  di¬ 
rection,  repeating  at  a  distance  equal  to  twice  the  primary  wavelength.  The  situation  of  detuned 
instability  is  between  the  fundamental  and  the  subharmonic  one. 

In  the  parabolized  stability  equations  (PSE)  approach  (Herbert  h  Bertolotti,  1987),  one  at¬ 
tempts  to  construct  an  approximate  solution  of  the  full  Navier-Stokes  equations,  where  the  nonpar¬ 
allel  effects  are  taken  into  account.  Though  the  parabolized  stability  equation  method  can  be  served 
as  an  alternative  for  DNS  in  some  range,  the  cost  of  PSE  increase  much  faster  than  that  of  DNS 
as  the  number  of  spanwise  modes  increases.  Therefore,  PSE  is  limited  in  practice  to  investigations 
of  narrow-band  focing  (Pruett  et  al.^  1995a). 

The  compressible  linear  stability  theory  (LST),  the  secondary  instability  theory  (SIT),  and 
the  parabolized  stability  equation  (PSE)  methods  can  only  provide  limited  information  for  flow 
transition,  and  they  are  far  away  from  practical  application. 

Direct  numerical  simulation  is  a  technique  that  all  the  flow  scales  are  accurately  resolved.  DNS 
of  flat-plate  boundary  layers  over  a  wide  range  of  Mach  numbers  (e.g.,  up  to  Mach  8  or  higher) 
have  resulted  in  some  encouraging  quantitative  comparisons  accurate  to  several  digits  with  theories 
(Pruett  et  a/.,  1995a).  However,  being  limited  by  the  capacity  of  modern  computers,  most  DNS  are 
limited  to  temporal  simulation  of  transition,  in  which  a  spatially  periodic  computational  domain 
travels  with  the  disturbance  and  the  temporal  evolution  of  the  disturbance  is  computed.  This 
enabled  simulations  into  the  later  stages  of  transition  (Zang  &  Hussaini,  1990;  Laurien  &  Kleiser, 
1989;  Pruett  h  Zang,  1992).  A  temporal  subharmonic  transition  in  a  Mach  4.5  flat-plate  boundary 
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layer  has  been  simulated  successfully  by  Adams  &  Kleiser  (1996)  using  direct  numerical  simulation. 
An  approach  has  been  developed  by  Liu  et  a/,  (1996a,  1996b)  to  simulate  the  whole  process  of 
transition  in  the  incompressible  boundary  layer  of  flat-plate  and  of  airfoils.  For  compressible  flow, 
some  of  the  techniques  from  their  previous  works,  are  combined  with  the  explicit  schemes,  which 
is  eflficient  especially  for  supersonic  flow  with  high  Mach  numbers.  In  addition,  the  explicit  code 
is  easy  to  be  vectorized  and  parallelized.  With  this  approach,  the  spatial  transition  of  boundary 
layer  flow  of  flat-plate  and  airfoil  is  investigated  by  direct  numerical  simulation  (Zhao  et  a/,  1997). 
However,  the  computational  cost  is  high  and  the  skin-friction  coefficients  and  velocity  profile  are 
not  accurate  for  fully  developed  turbulent  flow. 

Large  eddy  simulation  (LES)  is  a  technique  that  has  been  successfully  applied  to  the  study  of 
turbulent  flows,  and  has  become  a  very  important  method  of  simulation.  In  LES  only  the  large 
energy-carrying  scales  are  resolved,  while  the  influence  of  small  or  subgrid  scales  must  be  modeled 
appropriately.  A  filtering  process  is  usually  used  to  separate  the  large-  and  small-scale  motions, 
large-scale  structures  are  resolved  by  the  filtered  equations,  but  a  model  is  employed  to  formulate 
the  contributions  from  the  subgrid-scale  fluctuations,  such  as  the  subgrid  scale  stress  and  heat  flux 
terms.  Compared  with  DNS,  LES  is  more  affordable  at  present  stage.  LES  is  capable  to  solve 
problems  with  complex  geometry  using  spatial  approach  (Berlin,  1994;  Ducros  et  a/.,  1996). 

In  the  present  work,  the  spatial  transition  of  a  Mach  4.5  flat-plate  flow  is  investigated  by 
LES.  A  pair  of  oblique  first-mode  perturbation  is  imposed  on  the  inflow  boundary.  The  numerical 
method  is  based  on  the  three-dimensional  time-dependent  compressible  Navier-Stokes  equations 
in  the  curvilinear  coordinate  system.  A  compact  sixth-order  central  difference  scheme(Lele,  1992) 
is  applied  to  the  wall-normal  direction  and  streamwise  direction.  In  the  spanwise  direction,  the 
pseudo-spectral  method  is  used  based  on  the  assumption  that  the  periodic  condition  is  satisfied. 
The  compact  storage  third  order  Runge-Kutta  scheme  (Wray,  1986)  is  applied  for  time-integration. 
The  similar  approach  has  been  used  by  us  (Zhao  et  a/,  1997)  in  direct  numerical  simulation.  The 
major  objective  of  this  work  focus  on  the  LES  of  supersonic  flow.  The  numerical  results  that  are 
comparable  to  the  DNS  are  obtained  by  LES  with  less  grid  numbers  and  less  CPU-hours.  In  the 
present  work,  the  filtered  structure  function  model  (Ducros  et  al.^  1996)  is  adopted. 

The  paper  is  organized  as  follows.  Section  2  contains  the  governing  equations  and  the  LES 
model.  In  Section  3,  some  of  the  computational  details  are  presented.  The  computational  results 
and  discussion  are  in  Section  4. 


1.2  Mathematic  Model 

1.2.1  Governing  Equations 

In  large  eddy  simulation  of  fluid  flow  the  large  scale  structures  are  simulated  exactly  while  the 
effect  of  the  small  scale  structures  is  modeled.  The  large  scales  are  extracted  from  the  dependent 
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variables  by  applying  a  filtering  operation  to  the  continuity,  the  Navier-Stokes,  and  the  energy 
equations.  To  account  for  large  density  fluctuations  in  high-speed  compressible  flows,  the  Favre- 
filtering  operation  is  employed,  where  the  resolved  velocity  and  temperature  fields  are  written  in 
terms  of  Favre-filtered  quantities,  which  are  defined  as 

P 

where  the  ” — ”  denotes  the  spatial  filtering.  The  nondimensional  Favre-filtered  governing  equations 
of  continuity,  momentum,  and  temperature  are  described  as  follows: 


dpuk  ,9 _  dp  I  duki  ,  dTki 


(1.3) 


dpf  .  d 


dt  dx}z 


duk  ,  7(7 -l)M^  duk  d  'IP  df  dqk  ,  . 

^ki—  +  — ( a,.; )  +  ^  (1-4) 


Re 


dx{  dxk  Pr Re  dxk  dxk 


where  p  is  the  density,  Uk  is  the  velocity  component  in  the  kth  direction,  p  is  the  pressure,  and 
T  is  the  temperature.  The  viscous  stress  is 


^ki  =  P 


fduk  ^  dui')^  _  2dum 


\dxi  dxk)  3  5x 


4/ 


(1.5) 


In  the  nondimensionalization,  the  reference  values  for  length,  density,  velocity,  and  temperature  are 
^ini  Pooi  U-oo,  and  Too,  respectively.  is  the  displacement  thickness  of  inflow.  The  Mach  number, 
the  Reynolds  number,  the  Prandtl  number,  and  the  ratio  of  specific  heats,  are  defined  respectively 
as  follows: 


Moo  =  Uoo{lRToor^'\  Re 


Poo^oo^v 


Pr  = 


_  CpPo 


Cp 


poo  ^00 

where  R  is  the  ideal  gas  constant,  Cp  and  C'„  are  the  specific  heats  at  constant  pressure  and  constant 
volume.  Through  this  work,  Pr  =  0.7,  and  7  =  1.4.  The  viscosity  is  determined  according  to 
Sutherland’s  law,  in  dimensionless  form 


P  = 


+  5) 


S  = 


110.3/!: 


T-t-5  ’  - 

In  Eq.(1.3)  and  (1.4),  the  subgrid  scale  stress  and  heat  flux  are  denoted  by 

Tki  =  -P  (wfcu/  -  UkUl) 

Qk  =  -p  (ukT  -  UkT^ 


(1.6) 

(1.7) 


which  are  needed  to  be  modeled. 
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1.2.2  Filtered  Structure  Fuction  Model 


The  filtered  structure-function  model  is  developed  by  Ducros  et  al.  (1996).  The  subgrid  scale  shear 
stress  and  heat  flux  can  be  modeled  as 


Tkl  =  Pi't 


duic  dui\ 
dxi  dxk) 


2  dUjn 
ZdXm 


(1.8) 


yput  df 
Prt  dxk 


(1.9) 


Here,  Pr<  is  the  turbulent  Prandtl  number  taken  equal  to  0.6  as  in  isotropic  turbulence,  i/(  is  the 
turbulent  kinetic  viscosity  defined  as 


vtix,t)  =  0.0014C//^A  [Ff^(x,f)]^^^ 


(1.10) 


"  (3) 

where  F2  is  the  filtered  structure  function.  In  the  case  of  flat-plate  boundary  layer  flows  with 
meshes  flattened  in  the  wall-normal  direction,  F)  ’  takes  the  four-neighbor  formulation  proposed 
by  Normand  &  Lesieur  (1992). 


5(3) 


+ 


^t,3 


+ 


+ 


^»J-1 


(1.11) 


where 

A  =  is  used  to  characterize  the  grid  size.  Ck  is  the  Kolmogorov  constant  taking  the 

value  of  1.4.  is  a  discrete  Laplacian  filter  iterated  3  times,  which  is  served  as  a  high-pass 

filter  before  computing  the  structure  function.  The  first  iteration  of  the  Laplacian  filter  is 

defined  by 


(1.12) 


1.2.3  Boundary  Conditions 

Non-slip  boundary  condition  for  velocity  components  is  imposed  on  the  wall,  where  an  adiabatic 
condition  is  also  sited.  The  inflow  condition  is  composed  of  the  two-dimensional  Blasius-like  profile 
resulting  from  the  resolution  of  the  similarity  equation  and  a  pair  of  oblique  first-modes,  as  it  will  be 
described  in  detail  in  next  section.  At  the  far-field  boundary,  we  adopt  the  nonreflecting  boundary 
condition  proposed  by  Thompson(1987).  A  sponge  section  occupying  2  primary  wavelengths  in  the 
streamwise  direction  is  settled  near  the  outflow  boundary,  which  is  similar  to  the  sponge  condition 
used  by  Collis  &  Lele(1996).  A  source  term  is  added  to  the  right  hand  side  of  equations  as 

W{Q)^^f4x){Q) 
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where  the  sponge  function  is  given  by 


fd{x) 


X  e  (Xi.Xo] 
otherwise 


where  Xs  and  Xo  stand  for  the  streamwise  coordinate  of  the  starting  and  ending  points  of  the 
sponge  section,  respectively.  74^  and  Ns  are  used  to  control  the  amplitude  and  strength  of  the 
sponge  function.  Here,  ^4^  =  10  and  =  3  are  set. 


1.3  Computational  Procedure 

The  numerical  simulation  is  performed  using  a  spatial  approach  to  solve  a  full  compressible  Navier- 
Stokes  system  in  the  curvilinear  coordinates.  A  compact  sixth-order  central  difference  scheme 
(Lele,  1992)  is  applied  to  the  wall-normal  direction  and  streamwise  direction,  the  pseudo-spectral 
method  is  used  in  the  spanwise  direction.  The  compact  storage  third  order  Runge-Kutta  scheme 
(Wray,  1986)  is  adopted  for  time-integration.  The  computational  domain  is  displayed  in  Figure 
1.1.  We  designate  x  as  the  streamwise  direction,  y  as  the  spanwise  direction,  and  2;  as  the  wall- 
normal  direction.  The  flow  parameters  and  computational  details  have  been  collected  in  Table  1.1, 
where  the  inflow  displacement  thickness  is  denoted  by  Sin^  Xin  is  defined  as  the  distance  between 
the  leading  edge  of  the  flat-plate  and  the  upstream  boundary  of  the  computational  domain.  The 
Reynolds  number  is  based  on  the  displacement  thickness  at  the  inflow  boundary  and  the  free-stream 
velocity. 

The  computation  domain  covers  a  streamwise  region  between  x  =  91.335tn  and  x  =  510.2lJtn) 
thus  a  streamwise  dimension  of  =  418.886\n,  which  equals  to  32  streamwise  primary  wavelength 
Xx.  The  spanwise  dimension  is  Ly  =  T.blSin  covering  1  spanwise  primary  wavelength  Xy,  The 
wall-normal  dimension  at  the  upstream  boundary  is  228in.  The  number  of  grid  points  in  each 
direction  is  Nx  =  1024,  Ny  =  33,  and  Nz  =  61,  respectively.  The  grid,  uniformly  distributed  in  the 
streamwise  and  spanwise  direction,  is  stretched  from  the  wall  in  the  wall-normal  direction. 

A  pair  of  the  most  unstable  oblique  first  modes  with  wave-number  of  (1,1)  and  (1,-1)  ^nd  equal 
amplitude  are  imposed  as  the  inflow  perturbation.  The  frequency  of  the  disturbance  is  u;  =  0.3878. 
The  streamwise  and  spanwise  wavenumbers  are  =  0.48  and  j3  =  0.83  respectively,  corresponding 
to  the  streamwise  and  spanwise  wavelength  of  13.09  and  7.57,  thus  the  oblique  angle  is  60°.  The 
maxima  of  the  streamwise  velocity  disturbance  amplitude  is  about  2%  of  t/oo* 

The  time  step  is  set  to  be  0.0162  (1000  time  steps  are  used  per  period  of  the  forced  perturbation). 
With  the  number  of  grid  points  listed  in  Table  1.1,  the  CPU  time  of  each  time  step  is  8.2  seconds 
running  on  a  single  processor  of  Cray-T90.  For  this  simulation  we  used  30000  time  steps  which 
amounts  to  70  CPU  hours. 
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Figure  LI:  Computational  domain  of  flat-plate  boundary  layer  flow 

1.4  Results  and  Discussion 

1.4.1  Mean  flow  behavior 

In  this  section,  the  characteristics  of  the  mean  flow  is  discussed.  The  Fourier  transformation  is 
carried  out  in  time-  and  spanwise-direction  to  obtain  spectral  components  {kt^  ky)  in  Fourier  space, 
where  kt  and  ky  stand  for  the  temporal  and  spanwise  wavenumbers,  respectively.  The  mean  flow 
is  characterized  by  the  spectral  component  (0,0). 

The  skin  friction  coefficient  calculated  from  the  time-  and  spanwise-averaged  velocity  profile 
is  displayed  in  Figure  1.2.  The  spatial  evolution  of  skin  friction  coefficient  of  laminar  flow  is  also 
plotted  out  for  comparison.  It  is  observed  from  this  figure  that  the  sharp  growth  of  the  skin-friction 
coefficient  occurs  after  x  «  1805, n  (he.  Re^  —  1.8  X  10^),  which  will  be  defined  as  the  ’transition 
point’.  The  skin  friction  coeflScient  after  transition  is  in  good  agreement  with  the  flat-plate  theory 
of  turbulent  boundary  layer  by  Van  Driest(1956). 

The  thickness  of  boundary  layer  obtained  from  the  simulation  has  been  plotted  in  Figure  1.3. 
The  boundary  layer  thickness  are  calculated  based  on  the  time-  and  spanwise-averaged  velocity 
profile.  The  boundary  layer  thickness  is  defined  as  the  distance  from  the  wall  where  the  mean 
velocity  reaches  99%  of  the  free-stream  velocity.  A  sudden  growth  in  the  boundary  layer  thickness 
after  x  2005iri  Is  observed  as  a  result  of  the  transition  to  turbulence.  No  sudden  adjustment  can 
be  observed  from  the  curve  of  the  displacement  thickness  and  momentum  thickness. 

Time-  and  spanwise-averaged  streamwise  velocity  profiles  for  various  fixed  values  of  the  stream- 
wise  coordinate  x  =  91.335,n,  221. 635, n,  243. 385, 265.135i>,  and  286. 885, are  shown  in  Figure  1.4. 
The  inflow  .velocity  profile  at  a:  =  91.335,>,  is  of  a  typical  laminar  flow.  Starting  from  x  =  243. 385^^, 
the  mean  velocity  profiles  resemble  that  of  turbulent  flow. 
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Table  1.1:  Flow  parameters  and  computational  details 

In  Figure  1.5,  we  show  the  time-  and  spanwise-averaged  velocity  profile,  plotted  in  terms  of 
logarithm  scaled  wall  unit,  where  the  friction  velocity  and  the  friction  length  are  defined  as  Ur  = 
and  Zr  =  >  respectively.  The  linear  law  near  the  wall  has  been  marked  by  the 

cross  symbol.  The  log  law  curves  of  both  incompressible  and  compressible  turbulent  flow  are  also 
plotted  as  comparison.  The  wall  law  of  compressible  turbulent  flow  displayed  here  is  based  on 
the  effective  velocity  concept  of  van  Driest  (1951)  and  the  first-order  theoretical  law  of  wall  given 
by  White  and  Christoph  (1972).  Actually  the  curves  of  two  theories  coincide  very  well,  one  can 
hardly  distinguish  them  from  each  other  on  this  figure.  There  is  evidently  a  log  region  on  the 
profile  at  a:  =  243. 38(5, from  where  the  mean  velocity  profile  resembles  that  of  fully  developed 
turbulent  flow,  see  also  Figure  1.4.  In  the  log  region,  our  results  are  in  good  agreement  with  the 
theory  of  Van  Driest(1951)  and  White  &:  Christoph (1972).  The  lower  of  the  velocity  profile,  in  the 
log  range,  compared  with  that  of  incompressible  turbulent  flow,  can  be  explained  as  the  effect  of 
compressibility. 

1.4.2  Linear  and  nonlinear  disturbance  evolution 

In  Figure  1.6,  we  show  the  maxima  amplitude  of  streamwise  velocity  perturbation  of  selected  modes 
as  a  function  of  streamwise  coordinate.  The  wavenumbers  of  these  modes  are  denoted  as  {kt,ky), 
where  kt  refers  to  the  temporal  wavenumber,  ky  is  the  spanwise  wavenumber.  Here  and  after,  the 
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Figure  L2:  Streamwise  evolution  of  skin-  Figure  1.3:  Thickness  of  boundary  layer  as  a 

friction  coefficient  obtained  from  the  time-  function  of  streamwise  coordinate 

and  spanwise-averaged  streamwise  velocity 
profile 

modes  with  an  even  sum  of  kt  +  ky  are  called  the  even  modes^  those  with  an  odd  sum  of  kf  +  ky  are 
termed  as  the  odd  modes.  In  Figure  1.6,  the  even  modes  are  plotted  with  solid  line,  and  the  odd 
modes  with  dashed  lines.  It  is  apparent  from  the  figure  that  the  disturbance  growth  is  dominated 
by  even  modes. 

I 

As  can  be  seen  from  Figure  1.6,  at  the  inflow  boundary  (x*  =  91.33(J,n)j  the  superimposed 
primary  wave  (l,l)-mode  is  visibly  well  standing  out  above  all  other  modes.  Before  x  »  170^, 
the  (l,l)-mode  grows  with  its  linear  growth  rate,  which  is  approximately  a,*  =  -1.35  X  10“^. 
The  linear  growth  of  the  primary  wave  covers  almost  6  primary  wavelength  in  the  streamwise  di¬ 
rection.  It  is  noticeable  that  the  (0,2)-mode  is  rapidly  amplified  and  it  overtakes  the  primary 
wave  mode  at  x  150^in-  Since  the  primary  wave  is  composed  of  a  pair  of  oblique  modes 
(1,1)  and  (1,-1),  the  weakly  nonlinear  interaction  of  the  oblique  couple  can  excite  the  (0,2)and 
(2,0)-modes.  The  further  interaction  of  (l,l),(l,-l),{0,2),(2,0)-modes  results  in  the  generation 
of  (l,3),(3,l),(2,2),(3,3),(0,4),(4,0),(2,4),(4,2)-modes,  etc.  Thus,  the  weakly  nonlinear  interaction 
among  the  even  modes  can  only  excite  the  even  modes.  The  sharp  growth  of  the  odd  modes, 
which  are  evoked  by  the  strong  nonlinear  interaction,  is  observed  after  x  «  1505in-  If  Is  also  very 
interesting  to  see  that  the  linear  growth  of  primary  wave  extends  into  the  nonlinear  region.  At  far 
downstream,  both  the  even  modes  and  the  odd  modes  are  saturated,  as  the  maxima  amplitude  of 
velocity  perturbation  stays  near  a  certain  level. 

The  linear  and  nonlinear  disturbance  growth  before  x  «  ISO^tn  can  be  regarded  as  the  first 
stage  of  transition.  Up  to  this  point,  because  the  amplitude  of  disturbance  is  too  small  to  make 
evident  change  to  the  mean  flow  velocity  profile,  the  skin-friction  coefficient  stays  at  its  laminar 
flow  value.  The  study  of  the  growth  of  separated  Fourier  spectral  modes  may  provide  some  insights 
into  the  description  of  organized  motion  during  the  transition  process,  as  it  will  be  described  in 
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Figure  1.4:  Time-  and  spanwise-averaged  Figure  1.5:  Log-linear  plots  of  the  time- 

streamwise  velocity  profiles  at  different  and  spanwise-averaged  velocity  profile  in  wall 

streamwise  locations  unit.  The  log  law  of  incompressible  flat-plate 

boundary  layer  flow,  and  the  log  law  given 
by  Van  Driest(1951),  and  a  first-order  theo¬ 
retical  law  of  the  wall  given  by  White  and 
Chrisoph(1972)  are  also  plotted  for  compari¬ 
son. 

detail  in  the  next  section. 

1.4.3  Organized  motions 
Identification  of  A-vortices 

During  the  post-processing  of  our  LES  database,  there  are  several  possible  ways  to  identify  the 
large-scale  vortical  structures.  The  alternative  choices  include  the  low  pressure,  the  negative  value 
of  the  second  invariant  of  velocity  gradient  tensor,  and  the  peak  value  of  vorticity.  It  is  reported  by 
Sandham  &  Kleiser  (1992)  that  second  invariant  appears  to  be  a  good  measure  when  the  vortices 
are  very  weak,  and  the  low  pressure  can  get  better  results  for  strong  vortical  structures.  The 
vorticity  has  been  proved  to  be  problematic,  since  it  can  not  distinguish  between  shear  layers  and 
vortices.  It  is  also  shown  by  Kleiser  Laurien  (1985)  that  the  A-vortex  does  not  necessarily  consist 
of  a  buncL  of  vortex  line.  The  results  of  Robinson  (1991)  demonstrated  that  vortex  lines  can  give 
the  illusion  of  hairpin  vortices,  even  when  there  are  no  such  vortices  in  a  flow.  Thus  low  pressure 
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Figure  1.6:  Spatial  growth  of  the  maxima  amplitude  of  streamwise  velocity  perturbation  modes, 
the  even  modes  are  drawn  in  solid  line,  and  the  odd  modes  in  dashed  line. 

turns  out  to  be  the  most  useful  for  the  identification  of  vortices  than  other  criteria,  especially  at 
highly  nonlinear  stages(Adams  &  Kleiser,  1996).  In  our  work,  the  large-scale  vortices  are  identified 
by  the  iso-surface  of  low  pressure.  The  iso-surface  of  streamwise  vorticity  component  Ux  is  also 
used  as  a  creteria  for  comparison.  We  find  out  that  for  large-scale  vortical  structures,  low-pressure 
tube  usually  appeared  in  the  core  of  the  vortex.  The  large-scale  vortical  structures  can  be  identified 
precisely  through  both  the  low  pressure  criteria  and  the  concentration  of  streamwise  vorticity.  But 
the  results  obtained  from  these  two  methods  diversify  when  weakly  disturbed  flow  is  concerned, 
where  there  is  no  actually  physical  vortical  structure.  In  this  case,  although  the  iso-surface  of 
streamwise  vorticity  component  comes  out  to  be  tube-shaped,  they  can’t  be  interpreted  as  vortical 
structures.  For  example,  in  our  spatial  simulation,  an  oblique  wave  was  imposed  at  the  inlet 
boundary.  Near  the  inflow  boundary,  the  disturbance  experiences  a  linear  and  weakly  nonlinear 
growth.  In  the  region  between  x  «  91(5in  and  x  »  fhe  iso-surface  of  the  streamwise  vorticity 

component  appears  to  be  tube-shaped,  shown  in  Figure  1.7.  But  the  low  pressure  area  indicated 
by  iso-surface  of  pressure  in  this  figure  is  much  different  from  the  vorticity  tubes.  The  conclusion 
is  that  the  streamwise  vorticity  component  originating  from  the  initially  imposed  perturbation  can 
not  be  regarded  as  the  physical  vortical  structures.  It  is  also  quite  clear  from  this  figure  that  the 
axis  of  these  u>x  iso-surface  is  parallel  to  the  streamwise  direction.  So  the  low-pressure  method  is 
by  and  large  superior  compared  with  the  vorticity  method.  Hence,  the  low  pressure  will  be  used  to 
identify  the  vortical  structure. 

A  further  confirmation  of  the  low-pressure  method  in  our  results  is  depicted  in  Figure  1.8,  where 
the  iso-surface  of  low-pressure  (p  =  0.032pooU^)  is  plotted.  The  low-pressure  tubes  are  cut  by  two 
cross-section  planes.  On  each  plane,  the  streamline  of  cross-flow  is  drawn  as  well  as  the  contours  of 
streamwise  vorticity.  The  vortex  core  marked  by  the  streamline  just  sites  inside  the  low-pressure 
tube.  In  this  figure,  the  concentration  of  vorticity  coincides  with  the  low-pressure  part. 
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Figure  1.7:  Iso-surface  of  stream  wise  vortic- 
ity  and  low-pressure  near  the  inflow  bound¬ 
ary  {ux  =  ±0.06I7oo/5,„,  dark;  p  — 

0.035/)oot^,^  .light) 
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Figure  1.8:  Cross  flow  streamline  on  a  slice 
normal  to  the  streamwise  direction.  Iso¬ 
surface  (wire  form)  of  low  pressure  p  — 
O.Q32pooU^-  Contour  of  streamwise  vortic- 
ity  plotted  on  the  planes  of  cross-section 


Evolution  of  A- vortices 

The  Y-shaped  shear  layers  are  observed  in  temporal  simulation  of  subharmonic  breakdown  (Adams 
&  Kleiser,  1996)  and  oblique  breakdown  (Guo  et  al^  1994).  In  our  results  of  LES  using  the  spatial 
approach^  the  Y-shaped  shear  layer  being  displayed  by  the  iso-surface  of  spanwise  vorticity,  clearly 
appears  after  x  ^  170^tn*  The  spatial  distribution  of  these  Y-shaped  shear  layers  and  the  A- vortices 
is  depicted  by  a  top  view  in  Figure  1.9.  The  Y-shaped  shear  layers  are  staggered  in  the  streamwise 
direction.  On  the  cross-section  it  is  observed  that  the  rolling-up  of  the  Y-shaped  shear  layers 
generates  the  open-tip  A- vortices.  Looking  from  the  upstream,  the  right  branch  of  a  Y-shaped 
shear  layer  rolls  down  while  the  left  branch  of  the  staggered  Y-shaped  shear  layer  downstream  rolls 
up,  thus  a  clockwise  rotation  of  fluid  is  formed,  i.e.,  a  right  tail  of  an  open-tip  A-vortex.  Similarly, 
the  left  tail  of  the  open-tip  A-vortex  comes  from  the  rolling  down  of  a  left  Y-branch  and  rolling 
up  of  a  downstream  staggered  right  Y-branch,  Thus  the  resulted  open-tip  A-vortices  also  form  a 
staggered  system  in  the  streamwise  direction.  Because  of  the  deformation  and  break-up,  the  stalk 
of  the  Y-shaped  shear  layers  disappear  after  x  «  207 Jin,  while  the  remnant  of  the  Y-branches 
evolves  into  the  shear  layers  surrounding  the  tail  of  the  A-vortices. 

The  appearance  of  A-vortices  indicates  the  start  of  the  second  stage  of  transition.  A  series 
of  open-tip  A-vortices  between  a;  =  200Jin  a,nd  x  =  228Jin  a,re  displayed  in  Figure  1.10.  These 
A-vortices  are  staggered  in  the  streamwise  direction.  The  staggered  system  of  these  A-vortices 
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Figure  1.9:  Iso-surface  of  low  pressure  p  ==  0.0315/>oot^<^)  and  iso-surface  of  spanwise  vorticity 
Uy  =  0.68t/oo/<^m-  (data  duplicated  periodically  in  the  spanwise  direction). 

resembles  the  scenery  of  subharmonic  secondary  instability  (H-type)  in  a  flat-plate  boundary  layer 
flow,  the  H-type  breakdown  is  characterized  by  the  appearance  of  the  staggered  A-vortices.,  In 
the  case  of  a  subharmonic  secondary  instability,  the  inflow  perturbation  is  composed  of  a  two- 
dimensional  primary  wave  of  second-mode  and  random  noise.  The  wavenumber  {kfyky)  of  the 
primary  mode  is  (0,2).  Then  the  subharmonic  mode  with  a  wavenumber  of  (1,1)  grows  from 
the  random  noise.  The  strongly  nonlinear  interactions  between  the  dominant  mode  lead  to  the 
generation  of  a  staggered  system  of  A-vortices,  whose  streamwise  repeating  distance  equals  to 
twice  the  primary  wavelength.  Thus  streamwise  repeating  distance  is  the  same  as  the  wavelength 
of  the  subharmonic  mode. 

In  our  simulation,  a  pair  of  oblique  first-modes  with  wavenumber  (kt^ky)  of  (1, 1),  (1,  -1)  and 
equal  amplitude  are  imposed  on  the  upstream  boundary  to  serve  as  the  inflow  perturbation.  The 
instability  is  evoked  first  by  the  linear  and  weakly  nonlinear  interactions  between  the  oblique  modes 
and  the  even  modes,  then  by  the  strongly  nonlinear  interactions,  as  a  result  the  odd  modes  are 
brought  out.  As  the  onset  of  the  strongly  nonlinear  interactions  between  the  dominant  modes, 
the  staggered  A-vortices  appear.  The  streamwise  repeating  distance  of  the  A-vortices  is  close  to 
the  oblique  primary  wavelength.  Here,  the  wavenumber  of  the  primary  wave  is  comparable  to  the 
wavenumber  of  the  subharmonic  mode  in  H-type  instability.  Therefore,  the  A  vortical  structure 
coming  out  from  the  subharmonic  secondary  instability  and  from  the  oblique  instability  bear  some 
similarity,  as  it  is  described  by  Adams  &  Kleiser  (1993).  In  their  temporal  direct  numerical 
simulation  for  an  oblique  transition  at  Mach  number  4.5,  the  building-up  of  a  staggered  system 
of  streamwise  vortices  and  Y-shaped  shear  layers  is  observed  to  be  similar  to  the  subharmonic 
transition  (Guo  et  al,^  1994).  Although  the  mechanism  by  which  the  A-vortices  are  generated  is 
different  for  the  subharmonic  instability  and  oblique  instability,  the  resulted  A  vortical  structure  is 
quite  similar.  It  is  presumable  that  after  the  stage  of  A-vortices,  the  further  development  of  these 
two  types  of  transition  is  similar  to  each  other. 

Now  we  zoom  out  to  include  the  downstream  region  of  the  A-vortices  (see  Figure  1.11).  As  we 
have  seen  before,  from  x  =  2025in  to  x  =  2275,^,  the  tip  of  the  A-vortex  is  lifted  by  the  motion 
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Figure  1.10:  Iso-surface  of  low  pressure  p  =  0.03pooi^<^-  (data  duplicated  periodically  in  the 
span  wise  direction). 

induced  by  the  other  tail  of  the  same  vortex.  Although  the  lifting  angle  is  quite  small  (approximately 
15"^)  in  this  region,  the  downstream  velocity  difference  caused  by  the  mean  shear  drives  the  vortex 
to  stretch  considerably  in  the  streamwise  direction.  As  the  tip  of  A- vortex  sticks  out  of  the  strong 
mean  shear  layer,  the  decrease  in  the  wall-normal  derivative  of  the  streamwise  velocity  enhances 
the  inclination  of  the  tips  of  A-vortices.  (e.^.,  x  «  2345j’n,  x  ^  242^i„,  x  ^  2495in)*  The  A-vortices 
undergo  deformation  between  x  «  227 Sin  a-nd  x  ^  242Sin^  It  is  observed  from  the  side-view  of 
Figure  1.11  that  the  tips  of  the  A-vortex  lifts  dramatically  near  x  «  234(5in)  and  the  vortex  breaks 
into  two  parts  with  the  tip  and  tail  of  the  same  vortex  separating  each  other.  The  original  tail  of 
the  A-vortex  becomes  weak  as  the  the  iso-surface  of  low  pressure  shrinks.  The  broken  tails  can  be 
regarded  as  the  weak  vortex  torn  off  from  the  original  A-vortex,  but  they  are  still  visible  downstream 
near  x  «  24iSin>  The  connection  of  remnant  open-tips  of  the  original  A-vortex  makes  it  look  like 
the  head  of  a  hairpin  vortex  (near  x  ^  2‘i4Sin  on  Figure  1.11).  Near  x  «  240^tn  on  Figure  1.11, 
the  vortical  structure  is  more  like  hairpin  vortex  rather  than  the  A-vortex.  At  a:  «  24SSin,  the  A- 
vortices  have  eventually  evolved  into  the  hairpin  vortices.  The  head  of  the  hairpin  vortex  becomes 
strong,  as  the  iso-surface  of  low  pressure  grows.  The  legs  of  the  hairpin  vortex  are  now  still  broken 
with  the  head,  they  seem  to  be  the  remnant  of  the  broken  tails  of  the  A-vortices  at  x  230^,>i  and 
are  still  weak.  The  enclosed  angle  of  the  hairpin  vortex  head  is  about  53°.  The  appearance  of  the 
hairpin  vortex  is  regarded  as  the  third  stage  of  supersonic  boundary  layer  transition. 

Looking  from  the  upstream,  the  left  leg  of  a  hairpin-vortex  is  corresponding  to  negative  stream- 
wise  vorticity,  i.e.,  counterclockwise  rotation  {e.g.,  x  ^  2405{n^  y  ^  25,>;;  x  ^  248Sin^  y  ~  5(5, and 
X  ^  248(5,^,  y  ^  -25, on  Figure  1.11).  The  right  leg  with  a  positive  streamwise  vorticity  indicates 
a  clockwise  rotation (e. 5.,  x  2405^^,  z  -2.5in\  x  «  2485, n,  y  ~  25in;  and  x  ^  2485^,,,  y  «  -55,^ 
on  Figure  1.11).  The  small  vortices  staying  close  to  the  leg  of  hairpin  vortices  have  an  opposite 
sign  in  streamwise  vorticity. 
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Figure  1.11:  Iso-surface  of  low  pressure  p  =  O.O^pooU^,  (data  duplicated  periodically  in  the 
spanwise  direction). 

Both  vortices  and  shear  layers  are  interesting  flow  structures  of  the  transitional  and  turbulent 
flow.  The  A-vortices  and  the  evolving  hairpin  vortices  are  discribed  before  this  paragraph.  Now 
let  us  take  a  look  at  the  shear  layers  associated  with  the  vortical  structures.  In  Figure  1.12,  the 
contours  of  the  spanwise  vorticity  Uy  is  plotted  on  a  plane  cutting  through  the  vortex-tube  and  a 
symmetrical  plane  of  the  vortex  structure.  Here  the  spanwise  vorticity  Uy  is  a  good  measure  to  the 
shear.  The  shear  layers  locating  away  from  the  wall  are  produced  by  the  induced  motion  of  the 
vortical  structures. 

In  order  to  study  the  three-dimensional  relationship  between  the  shear  layers  and  the  vortical 
structures,  the  contours  of  spanwise  vorticity  is  plotted  on  several  cross-sections  through  x  «  200(5jn 
and  X  «  24bSin  in  Figure  1.13.  It  is  observed  that  each  tail  of  a  A-vortex  is  wrapped  by  two  high 
shear  regions,  i.e.,  one  top  layer  and  one  bottom  layer.  Both  the  upper  and  lower  shear  layers  are 
the  branches  of  the  original  Y-shaped  shear  layer  shown  before  in  Figure  1.9. 

The  interactions  between  the  shear  layers  and  the  vortical  structures  are  very  complex.  Gen¬ 
erally  speaking,  the  first  appearance  of  the  A-vortices  can  be  explained  as  a  result  of  shear  layer’s 
rolling  up.  And  new  shear  layers  are  generated  around  the  vortices  by  a  process  of  vortex  stretch¬ 
ing  and  convection  (Stuart  1965).  Usually  two  types  of  shear  layers  are  generated  by  the  quasi- 
streamwise  vortical  structure,  the  detached  and  the  attached  shear  layer.  Only  the  detached  shear 
layer  is  discribed  in  the  previous  paragraph,  the  attached  shear  layer  is  mixed  with  the  mean  shear. 
The  attached  high-shear  layer  is  visible  by  plotting  the  contours  of  the  wall  shear  stress,  as  seen  in 
Figure  1.14,  where  the  wall  shear  stress  distribution  on  the  wall  in  the  region  between  x  =  2005, n 
and  X  =  2505in  is  displayed.  High  wall  shear  rate  corresponding  to  the  attached  shear  layer  can  be 
found  close  to  the  outside  of  the  A-vortex  tail,  e.flf.,  x  «  2225, n,  y  ^  -5,  —3, 3, 55,^,  or  a;  «  2285,^, 
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Figure  1.12:  Iso-surface  of  low  pressure  p  =  O.OSpooU^,  contours  of  instantaneous  spanwise  vorticity 
on  the  streamwise-sections.  (data  duplicated  periodically  in  the  spanwise  direction). 


Figure  1.13:  Iso-surface  of  low  pressure  p  =  contours  of  instantaneous  spanwise  vorticity 

on  the  cross-sections,  (data  duplicated  periodically  in  the  spanwise  direction). 


y  ^  -7,  -1, 1, 76in-  The  attached  shear  layer  is  generated  by  the  down-wash  motion  of  fluid  toward 
the  wall  induced  by  the  vortex  tail.  The  opposite  position  of  these  high  shear  region  with  respect  to 
the  vortex  tail  is  controlled  by  the  up-wash  motion,  where  the  wall  shear  rate  is  low.  The  strongest 
wall  shear  is  discovered  between  the  two  tails  of  the  open-tip  A-vortices  or  the  two  legs  of  the 
hairpin  vortices. 
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Figure  1.14:  Iso-surface  of  low  pressure  p  =  0.03pooU^,  contours  of  the  wall  shear  stress,  (data 
duplicated  periodically  in  the  spanwise  direction). 

The  further  development  of  these  hairpin  vortices  is  displayed  in  Figure  1.15.  In  this  region,  the 
spanwise  expansion  of  the  hairpin  vortex  is  quite  evident,  since  the  low-pressure  areas  of  adjacent 
hairpin  vortices  tend  to  merge  together.  Actually  the  vortices  have  an  opposite  rotating  direction. 
A  deformation  of  the  hairpin  vortices  comes  out  to  be  visible  in  this  figure,  and  the  legs  of  the 
hairpin  vortices  are  stretched  in  the  streamwise  direction,  as  the  vortical  structure  continues  to 
expand  in  the  spanwise  direction.  The  spanwise  distance  between  the  left  and  right  legs  of  the 
vortex  at  a;  3205tn  reaches  6.45in*  The  legs  of  the  hairpin  vortices  are  almost  parallel  to  each 
other,  and  also  parallel  to  the  wall.  Besides  the  deformation,  these  hairpin  vortices  also  become 
weaker  as  the  low-pressure  areas  shrink.  The  head  of  the  vortical  structures  are  almost  invisible. 
The  broken  of  these  hairpin  vortices  may  lead  to  the  generation  of  some  small  scale  structures  near 
the  head  of  the  former  hairpin  vortices.  But  they  are  not  observed  clearly  in  our  simulation,  which 
may  be  blamed  to  lack  of  resolution  in  this  region.  Further  high-resolution  DNS  may  be  required. 
The  visible  structures  left  after  the  breakdown  are  the  streamwise  vortices  which  are  located  very 
close  to  the  wall.  They  evolve  from  the  legs  of  the  former  hairpin  vortical  structure. 


1.5  Conclusions 

The  transition  process  in  a  Mach  number  4.5  flat-plate  boundary  layer  is  investigated  by  large  eddy 
simulation  using  the  spatial  approach.  A  pair  of  oblique  first-mode  perturbation  is  imposed  at  the 
inflow  boundary,  and  the  initial  field  is  a  laminar  boundary  layer  flow. 
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Figure  1.15:  Iso-surface  of  low  pressure  p  =  O.OSpooU^.  (data  duplicated  periodically  in  the 
spanwise  direction). 

The  linear  and  weakly  nonlinear  growth  region  covers  from  x  «  9l4„  to  a;  «  1705, The  spatial 
growth  rate  of  the  primary  wave  remains  at  the  linear  value  until  x  «  1705i„.  Because  the  primary 
wave  is  composed  of  a  pair  of  oblique  modes  (1, 1)  and  (1,-1),  only  the  ’even’  modes(with  an  even 
sum  of  wavenumbers  kt  -f  ky)  can  be  excited  through  the  weakly  nonlinear  interaction. 

The  strongly  nonlinear  interaction  occurs  after  x  fa  1705i„,  which  is  characterized  by  the  ap¬ 
pearance  of  a  staggered  system  of  A-vortices  and  Y-shaped  shear  layers.  The  ’odd’  modes  (with  an 
odd  sum  of  wavenumbers  kt  +  ky)  are  evoked  by  the  strongly  nonlinear  interaction. 

Shear  layers  and  vortices  are  the  two  major  structures  in  the  transitional  and  turbulent  flow. 
Typical  A-vortices  are  visible  between  x  fa  1805, •„  and  x  fa  2265, The  vortex  core  of  the  A- 
vortices  are  wrapped  by  the  branches  of  the  Y-shaped  shear  layers.  The  induced  motion  of  each 
tail  working  on  the  other  gives  the  inclined  nature  of  the  vortices.  In  the  downstream  part  of  this 
region,  the  deformation  of  the  A-vortices  and  the  Y-shaped  shear  layers  are  observed.  As  a  result 
of  the  interaction  between  the  self  induced  motion  and  the  mean  shear,  deformation  occurs  to  the 
A-vortices  and  the  Y-shaped  shear  layers. 

The  next  stage  of  transition  is  characterized  by  the  appearance  of  the  hairpin  vortices,  which 
first  appear  at  a;  «  2485, The  hairpin  vortices  are  evolved  from  the  deformation  of  the  open- tip 
A-vortices.  The  further  growth  of  the  hairpin  vortices  is  observed  downstream,  where  the  spanwise 
scale  of  the  hairpin  is  extended  with  increased  inclination.  The  hairpin  vortices  are  matured  near 
X  «  3005,„  where  the  strength  of  the  vortex  reaches  a  maximum.  After  x  fa  3105,„,  the  hairpin 
vortices  experience  a  considerable  deformation,  the  legs  stretch  in  stream  wise  direction  and  evolve 
into  the  streamwise  vortices  parallel  to  the  wall.  The  heads  of  the  hairpin  vortices  become  weak, 
which  is  denoted  by  a  thin  or  discontinued  strand  of  low  pressure. 
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The  tails  of  the  A-vortices  and  the  legs  of  the  hairpin  vortices  generate  detached  shear  layers 
cis  well  as  high  wall  shear.  The  detached  shear  layers  roll  up  into  new  vortices. 

The  above  results  are  in  good  agreement  with  the  direct  numerical  simulation  both  with  a 
temporal  approach  (Adams  &  Kleiser,  1996)  and  a  spatial  approach  (Zhao  e/,  1997).  Compared 
with  our  previous  DNS  work  (Zhao  e/,  1997),  fewer  grid  nodes  and  computational  time  are  required 
by  LES,  but  the  transitional  process  is  simulated  accurately.  The  mean  flow  features  including  the 
skin  friction  coefficient,  log  law  distribution  of  mean  velocity  are  agreed  well  with  the  theoretical 
results.  Most  of  the  important  organized  structures  are  also  observed  in  the  LES  data  base,  which 
is  helpful  to  the  understanding  of  the  basic  mechanism  of  the  flow  transition.  This  approach  has 
also  been  applied  to  the  LES  of  supersonic  and  subsonic  flow  around  airfoils,  as  it  will  be  reported 
later. 


Chapter  2 


Compressible  Flow  Around  Two-Dimensional  Airfoil 


2.1  Two-Dimensional  Grid  Generation 

An  elliptic  grid  generation  method  first  proposed  by  Spekreuse  (1995)  is  used  to  generate  2D 
grids.  The  elliptic  grid  generation  method  is  based  on  a  composite  mapping,  which  is  consisted 
of  a  nonlinear  transfinite  algebraic  transformation  and  an  elliptic  transformation.  The  algebraic 
transformation  maps  the  computational  space  C  onto  a  parameter  space  P,  and  the  elliptic  trans¬ 
formation  maps  the  parameter  space  on  to  the  physical  domain  P.  The  computational  space, 
parameter  space,  and  the  physical  domain  are  illustrated  in  Figure  2.1. 

^  t 

E4  E4 

^ j - 

El  E2  El  E2 

ol _ , _ _  ol _ 

0 E3 t  0  E3  1  s 

Computational  Parameter 

Space  Space 
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Figure  2.1:  Computational  space  C,  Parameter  space  P,  and  Physical  domain  D 

The  computational  space  C  is  defined  as  the  unit  square  in  a  two-dimensional  space  with  Carte¬ 
sian  coordinates  (^,  rj),  and  ^  €  [0, 1],  r]  6  [0, 1]  (see  Figure  2.1).  The  grids  are  uniformly  distributed 
on  the  boundaries  and  in  the  interior  area  of  the  computational  space.  The  mesh  sizes  are 
in  the  ^  direction  and  in  the  i]  direction,  where  and  Nr,  are  the  grid  numbers  in  the 

corresponding  direction.  The  parameter  space  V  is  defined  as  a  unit  space  in  a  two-dimensional 
space  with  Cartesian  coordinate  (s,  t),  and  s  €  [0, 1],  t  G  [0, 1].  The  boundary  values  of  s  and  t  are 
determined  by  the  grid  point  distribution  in  the  physical  domain. 
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•  s  =  0  at  edge  Ei  and  s  =  1  at  edge  E2 

•  s  is  the  normalized  arclength  along  edges  £3  and  £4 

•  t  =  0  at  edge  £3  and  t  =  1  at  edge  £4 

•  t  is  the  normalized  arclength  along  edges  £1  and  £2 

An  algebraic  transformation  $  :  C  V  is  defined  to  map  the  computational  space  C  onto 
the  parameter  space  V.  The  grid  distribution  is  specified  by  this  algebraic  transformation,  which 
depends  on  the  prescribed  boundary  grid  point  distribution.  The  interior  grid  point  distribution 
inside  the  domain,  generated  by  the  algebraic  transformation,  is  a  good  reflection  of  the  prescribed 
boundary  grid  point  distribution.  Let  and  se^(^)  =  s(C)  1)  denote  the  normalized 

arclength  along  edges  £3  and  £4,  and  tEiiv)  =  denote  the  normalized 

arclength  along  edges  £1  and  £2.  The  algebraic  transformation  s  :  C  —)■  P  is  defined  as 


S  =  S£;3(0(1-0 

Equation  (2.1)  is  called  the  algebraic  straight  line  transformation.  It  defines  a  differentiable  one- 
to-one  mapping  because  of  the  positiveness  of  the  Jacobian:  s^trj  ~  Sjjt^  >  0. 

The  elliptic  transformation  x  :  V  which  is  independent  of  the  prescribed  boundary  grid 

point  distribution,  is  defined  to  map  the  parameter  space  V  onto  the  physical  domain  V,  The 
elliptic  transformation  is  equivalent  to  a  set  of  Laplace  equations 


^xx  d"  ^yy  — 

^xx  ^yy  ” 

The  elliptic  transformation  defined  by  the  above  equations  is  also  differentiable  and  one-to-one. 

Till  now  we  have  defined  two  transformations,  i.e.,  the  algebraic  transformation  s  :C  and 
the  elliptic  transformation  x  :V^V.  Because  both  the  algebraic  transformation  and  the  elliptic 
transformation  are  differentiable  and  one-to-one.  The  composition  the  two  transformation  is  also 
differentiable  and  one-to-one,  so  as  to  the  inverse  transformation. 

In  physical  domain,  the  curvilinear  coordinate  system  satisfies  a  system  of  Laplace  equations: 

Ar  0  (2.3) 

where  r  =  (x,y)^.  The  inherent  smoothness  of  the  Laplace  operator  makes  the  grids  smoothly 
distributed  in  the  physical  domain.  Being  transformed  to  the  computational  space,  this  Laplace 


0 

(2-2) 
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system  becomes  a  set  of  Poisson  equations.  The  control  functions  is  determined  by  the  composed 
transformation  according  to  the  following  procedures. 

First,  Eq.(2.2)  is  transformed  into  the  computational  space  C: 


As  =  +  2g^‘^s^r,  +  +  Agsr, 

At  =  +  2g^'^t^r)  +  g^^trjr,  +  A^t^  +  Agt^, 


where  are  the  components  of  the  contravariant  metric  tensor,  which  can  be  calculated 

from  the  covariant  metric  tensor 

</"  =  ^322  =  {Vr^rr,) 

=  -^512= 

J  is  the  Jacobian:s^t^  -  From  Eq.(2.2)  and  (2.4),  we  have 


where 


A^ 

Aij 


and  the  matrix  T  is  defined  as 

tn  ) 

Then  the  Laplace  system  Eq.(2.3)  is  transformed  to  the  computational  space  C: 

+  2^^  +  A^r^  +  AgVr,  =  0, 


(2.5) 


(2.6) 


(2.7) 


(2.8) 


Substitute  Eq.(2.5)  into  Eq.(2.8),  and  A?;  are  replaced  by  the  control  functions  on  the  right- 
hand-side  of  Eq.(2.5),  and  we  obtain  the  Poisson  equations  for  the  grid  generation  as  follows: 


where  the  control  functions  -Pn^ ^22^ ^22^  determined  by  the  algebraic  trans¬ 
formation,  as  defined  previously  in  Eq.(2.6). 

The  elliptic  transformation  is  carried  by  solving  a  set  of  Poisson  equations.  The  control  functions 
are  specified  by  the  algebraic  transformation  only  and  it  is,  therefore,  not  needed  to  compute  the 
control  functions  at  the  boundary  and  to  interpolate  them  into  the  interior  of  the  domain,  as 
required  in  the  case  for  all  well-known  elliptic  grid  generation  systems  based  on  Poisson  systems. 
The  computed  grids  are  in  general  not  orthogonal  at  the  boundary.  The  algebraic  transformation 
can  be  redefined  to  obtain  a  grid  which  is  orthogonal  at  the  boundary. 

We  developed  a  numerical  grid  generation  code  based  on  this  method.  Any  structured  grid  with 
four  boundary  edges  can  be  generated,  the  orthogonality  condition  can  also  be  imposed  at  any  edges. 
The  continuity  of  grid  intervals  in  the  interior  area  is  guaranteed  by  the  elliptic  equations,  while 
the  continuity  of  grid  intervals  on  the  boundary  edges  is  set  by  users.  In  order  to  achieve  the  best 
results,  the  boundary  grid  distribution  usually  should  be  smooth  (continuous  up  to  second  order 
of  derivative).  Here  we  provide  some  examples  of  2D  numerical  generation.  Figure  2.2  shows  the 
C-grid  around  a  Joukowsky  airfoil.  The  grids  near  the  leading  edge  and  trailing  edge  are  displayed 
in  Figure  2.3,  which  shows  that  the  smoothness  and  orthogonality  are  well  maintained.  The  grid 
number  is  741  in  streamwise  (^)  direction  and  121  in  wall-normal  (17)  direction. 


Figure  2.2:  C-grid  around  a  Joukowsky  airfoil. 

In  Figure  2.4,  the  C-grid  around  a  NACA  0012  airfoil  is  displayed.  The  grids  are  orthogonal 
at  the  boundaries.  The  distribution  of  the  grids  near  the  leading-edge  and  trailing-edge  are  also 
depicted  by  Figure  2.4(c)  and  (d),  where  good  orthogonality  has  also  been  achieved,  even  though 
there  are  two  singular  points  at  the  trailing-edge.  The  number  of  the  grids  is  601  in  the  streamwise 
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(a)  grids  near  the  leading  edge  (b)  grids  near  the  trailing  edge 

Figure  2.3:  C-grid  near  the  leading  edge  and  trailing  edge  of  a  Joukowsky  airfoil.  (The  grids  are 
orthogonal  on  the  airfoil  surface) 

(^)  direction  and  121  in  the  wall-normal  (tj)  direction. 
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(c)  grids  near  the  leading-edge 


(d)  grids  near  the  trailing  edge 


Figure  2.4:  C-grid  around  a  NACA  0012  airfoil 
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2.2  NACA  0012  Airfoil  with  a  20°  Angle  of  Attack 


The  situation  of  an  airfoil  with  a  high  angle  of  attack  occurs  in  many  aerodynamic  applications 
which  are  of  particular  interest.  In  these  cases,  the  fluid  flow  around  the  airfoil  becomes  very 
unstable  and  eddy  structures  are  formed  in  the  vicinity  of  the  airfoil.  These  eddy  structures  lead  to 
modification  of  the  airfoil  aerodynamic  coefficients  which  can  produce  noises.  In  order  to  reproduce 
the  unsteady  development  of  the  eddy  structures  in  the  vicinity  of  airfoils,  a  large  eddy  simulation 
has  been  launched  for  the  compressible  flow  around  a  NACA  0012  airfoil  at  a  20  degree  angle  of 
attack.  The  Reynolds  number  Re  =  5  x  10^  based  on  the  chord  length  and  freestream  velocity. 
The  Mach  number  is  Af  =  0.4.  A  C-grid  is  generated  using  the  method  introduced  in  the  previous 
section.  The  grid  number  is  601  in  the  ^-direction  and  121  in  the  77-direction. 

The  results  are  shown  in  Figure  2.5,  where  the  contours  of  spanwise  vorticity  at  different  time 
are  displayed.  In  the  first  frame  (Figure  2.5(a)),  a  separation  bubble  is  observed  near  the  leading 
edge.  The  leading-edge  separation  bubble  is  formed  when  the  laminar  boundary  layer  separates 
from  the  surface  as  a  result  of  the  strong  adverse  pressure  gradient  downstream  of  the  point  of 
minimum  pressure.  This  minimum-pressure  point  moves  upstream  as  the  angle  of  attack  increases, 
while  the  length  of  the  separation  bubble  decreases  (Arena  and  Mueller,  1980).  In  our  case,  the 
angle  of  attack  is  20  degree,  the  separation  bubble  near  the  leading-edge  is  very  short.  Because  the 
separation  bubble  is  very  unstable,  eddies  keep  shedding  from  the  leading-edge.  Three  large  scale 
eddy  structures  can  be  observed  in  Figure  2.5(a).  in  frame  (b)  of  Figure  2.5,  the  most  downstream 
eddy  reattaches  on  the  airfoil  surface.  A  new  shear  layer  is  formed  near  the  surface  as  a  result  of 
the  eddy  reattachment.  In  frame  (c)  the  reattached  eddy  merges  with  the  second  and  the  third 
eddy,  the  three  vortices  is  rolling  into  one  big  vortex.  The  induced  shear  layer  evolves  into  a  strong 
vortex  with  an  opposite  rolling  direction.  The  vortex-couple  with  opposite  rolling  direction  can 
be  observed  clearly  in  frame  (e)  of  Figure  2.5.  In  the  same  time,  new  eddy  structure  is  shedding 
from  the  leading-edge  separation  bubble.  The  vortex-couple  is  rolling  downstream  as  shown  in 
frame  (e)-(g),  until  it  hits  the  wake,  where  some  small  eddy  structures  exist.  In  Figure  2.6,  a 
zoom-in  area  near  the  trailing  edge  is  displayed,  each  frame  is  corresponding  to  the  same  frame  in 
Figure  2.5.  In  frame  (a)  of  Figure  2.6,  instability  is  observed  near  the  trailing  edge.  As  a  result 
of  wake  instability,  alternate  eddy  structures  are  generated  in  frame  (d)  and  (e).  The  interaction 
between  the  large  scale  eddy  structure,  e.g.  the  vortex-couple  and  the  small  scale  vortices  inside 
the  wake,  is  displayed  through  frame  (h)  to  (i)  of  Figure  2.6.  As  a  matter  of  fact,  the  large  scale 
eddy  structures  Overtake  the  small  scale  vortices,  which  disappeared  after  the  interaction,  and  new 
large  scale  vortex  is  induced. 


26 


Spanwse  Vorticity  (ci) 
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Figure  2.5:  Contours  of  spanwise  vorticity  at  different  stage 
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Figure  2.6:  Contours  of  spanwise  vorticity  near  the  trailing  edge  at  different  stage 
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Chapter  3 


Flow  Around  Three-Dimensional  Airfoil 


3.1  Three-Dimensional  Grid  Generation 

The  basic  idea  of  three-dimensional  grid  generation  is  similar  to  that  of  the  two-dimensional  case. 
The  computational  space  is  a  unit  cubic  with  ^  £  [0, 1],  rj  6  [0, 1],  C  €  [0, 1].  The  parameter  space  is 
a  unit  cubic  with  s  £  [0, 1],  t  £  [0, 1],  u  £  [0, 1],  see  Figure  3.1. 

•  s  =  0  at  face  Fi  and  s  =  1  at  face 

•  s  is  the  normalized  arclength  along  edges  £'i,  F2,  £'3  and  £4 

•  t  =  0  at  face  £3  and  t  =  1  at  face  £4 

•  t  is  the  normalized  arclength  along  edges  £5,  £6,  £7  and  £3 

•  u  =  0  at  face  £5  and  t  =  1  at  face  Fq 

•  t  is  the  normalized  arclength  along  edges  £9,  £10,  £11  and  £12 

Let  =  s(^,0,0),  SE^iO  =  SE^iO  =  s(^,0, 1),  and  se,{0  =  1. 1)  denote 

the  normalized  arclength  along  edges  £1,  £2,  £3,  and  £4;  tE^iri)  =  t(0, 7?, 0),  tEeii])  =  0), 

^Ejiv)  —  ^(0,7?,  1),  and  tE^iri)  =  1)  denote  the  normalized  arclength  along  edges  £5,  Ee,  £7 

and  £3;  UE^iO  =  u{0,QX),UEioiO  =  ^(l.O.C),  wEu(C)  =  “(O.  l,C),and  UEi^iQ  =  7i(l,  1, 77)  denote 
the  normalized  arclength  along  edges  £"9,  £10,  £11  £12.  The  algebraic  transformation  from 

computational  space  to  parameter  space  is  defined  as 


s  =  SEi  (0(1  -  0(1  -  ■‘i)  +  sf;2 (0^(1 -“) +  «jE3 (0(1  -  0w  + 

t  =  tEi{v){l- s){l-u)  +  tEeiv)s{'^-u)  +  tET{r])i'i-- s)u  +  tE^{ri)su  (3.1) 

u  =  UEg{0{l-s){l-t)  +  UE^oiOs{l-t)  +  UE^^{0{l-s)t  +  UE^^{Ost 
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Figure  3.1:  Computational  space  Parameter  space  P,  and  Physical  domain  D 


Equation  3.1  is  called  the  algebraic  straight  line  transformation. 

The  Poisson  equations  for  the  grid  generation  are  as  follows 

+  {g'^Ph  +  g^^Ph  +  g^^Pk  +  ^g'^Ph  +  ^g'^Ph  +  ^g^^PL)ri 

+  {g^^Ph  +  g^'^Ph  +  g^^Piz  +  ‘^g^^Ph  +  "^g'^Ph  +  ‘^g^^PhVr, 

+  (ff^‘Pn+5^^P|2  +  5^^^33  +  2<7*^^2  +  25^^if3  +  2f7^^if3)^C  =  0  (3-2) 

where  are  contravariant  metric  tensor,  which  are  calculated  through  the 

covariant  metric  tensor 

5“  =  ^(fl22<733  -  5I3) 

5^^  =  ^(5ii<733  -P?3) 

g^^  =  ^(5llff22  -  <7?2) 
g^^  =  j2  (^i3ff23  -  912933) 

g^^  =  -j2{9l2923  -  913922) 

5^^  =  ^{g\29i3-  923911) 


where  control  functions  are  defined  as 


(sa] 

^  Sr,v  ^ 

Pn  =  -T-^ 

ke 

,  P22  =  -T-^ 
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He  ) 

\  '^VV  j 
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kc 
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Pi2  =  -T-' 

f '^"1 

hn 

,  Pi3  =  -T-' 

(  s,c] 
kc 

,  P23  =  -T-' 

(Sr,c\ 

J 

) 

and  the  matrix  T  is  defined  as 


tf  tri  t/* 


V  H 


H  j 


(3.3) 


(3.4) 


Figure  3.2  shows  the  C-H-grid  around  a  3D  delta  wing,  where  the  grids  on  the  wing  surface  and 
on  a  ^  -  C  plane  are  displayed.  We  choose  C-type  for  the  streamwise  direction  and  H-type  for  the 
spanwise  direction. 


Figure  3.2:  C-H-grid  around  a  3D  delta  wing. 

A  set  of  H-C-grid  around  a  80°  delta  wing  is  generated  with  the  same  method.  In  this  case,  we 
choose  H-type  for  the  streamwise  direction,  and  C-type  for  the  spanwise  direction.  The  grids  are 
also  required  to  be  orthogonal  at  the  boundaries.  The  grids  on  selected  ^  —  f]  sections  are  displayed 
in  Figure  3.3(a).  The  grids  are  orthogonal  on  delta  wing  surface  and  boundaries.  Both  coarse  and 
fine  grid  are  generated,  and  the  fine  grid  (141  x  70  x  70)  is  shown  in  Figure  3.3(b). 


3.2  3D  Delta  Wing  with  a  20°  Angle  of  Attack 

Some  preliminary  results  of  the  large  eddy  simulation  of  compressible  flow  around  a  3D  delta  wing 
with  a  20°  angle  of  attack  are  obtained.  The  Reynolds  number  based  on  the  maximum  chord  length 
and  the  freestream  velocity  is  5  x  10®.  March  number  is  M  =  0.4.  The  grid  number  along  the  ^-,  rj- 
and  C-direction  are  221,  111  and  71,  respectively.  Here  the  streamwise,  spanwise,  and  wall-normal 
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(a)  coarse  grid 


(b)  fine  grid 
Figure  3.3:  H-C-grid  around  a  80®  delta  wing 

directions  are  denoted  by  77,  and  ^  respectively.  The  computational  domain  only  contains  half  of 
the  delta  wing  and  symmetric  condition  is  imposed  on  the  symmetric  plane. 

In  Figure  3.4,  the  contours  of  instantaneous  pressure  on  selected  ^  —  C  sections  are  displayed. 
The  low  pressure  areas  are  generally  corresponding  to  flow  separation  and  eddy  structures  shedding 
from  the  wing  surface. 


Figure  3.4:  Contours  of  instantaneous  pressure  on  selected  ^  -  C  planes. 

Three  dimensional  streamlines  starting  from  the  leading-edge  of  the  delta  wing  are  displayed  in 


Figure  3.5.  From  this  figure,  several  vortex  system  can  be  recognized.  A  strong  streamwise  vortex 
starts  from  the  middle  of  the  wing.  This  vortex  system  comes  from  the  separation  near  the  wing 
surface,  and  it  tends  to  move  outward  to  the  symmetric  plane.  The  wing-tip  vortices  can  also  be 
observed,  and  it  tends  to  move  toward  the  symmetric  plane. 


Figure  3.5:  3D  streamline  starting  from  the  leading-edge 

Three  components  of  the  vorticity  are  displayed  in  Figure  3.6.  Again,  the  separation  and 
shedding  of  eddy  structures  near  the  wing  surface  can  be  read  from  this  figure.  It  is  also  observed 
that  the  vortex  system  near  the  wing-tip  is  very  complex  and  it  may  require  more  resolution  than 
we  have  in  this  case.  The  strong  streamwise  vortices  observed  in  Figure  3.5  is  also  quite  clear  by 
looking  at  the  contours  of  streamwise  vorticity  in  Figure  3.6(a). 
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(b)  Contours  of  spanwise  vorticity 


(c)  Contours  of  wall-normal  vorticity 
Figure  3.6:  Contours  of  vorticity  on  selected  ^  (  planes. 
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Figure  2.5:  Contours  of  spanwise  vorticity  at  different  stage 


Figure  2.6:  Contours  of  spanwise  vorticity  near  the  trailing  edge  at  different  stage 


(a)  Contours  of  streamwise  vorticity 


(c)  Contours  of  wall-normal  vorticity 


Figure  3.6:  Contours  of  vorticity  on  selected  plane 


